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The PT symmetric potential Vo[cos(27rx/a) + iX sin(27rx/a)] has a completely real 
spectrum for A < 1, and begins to develop complex eigenvalues for A > 1. At the 
symmetry-breaking threshold A = 1 some of the eigenvectors become degenerate, giv- 
ing rise to a Jordan-block structure for each degenerate eigenvector. In general this 
{Sj . is expected to give rise to a secular growth in the amplitude of the wave. However, 
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it has been shown in a recent paper by Longhi, by numerical simulation and by the 
use of perturbation theory, that for an initial wave packet this growth is suppressed, 
£f) ' giving instead a constant maximum amplitude. We revisit this problem by devel- 

oping the perturbation theory further. We verify that the results found by Longhi 
persist to second order, and with different input wave packets we are able to see the 
seeds in perturbation theory of the phenomenon of birefringence first discovered by 

Oh" Makris et al. 

O ' 

c/3 1 PACS numbers: 42.25.Bs, 02.30.Gp, 11.30.Er, 42.82.Et 
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I. INTRODUCTION 



The study of quantum mechanical Hamiltonians that are PT-symmetric but not 
Hermitian[l|-j6j has recently found an unexpected application in classical optics[3]-[l5[, due 
to the fact that in the paraxial approximation the equation of propagation of an electromag- 
netic wave in a medium is formally identical to the Schrodinger equation, but with different 
interpretations for the symbols appearing therein. It turns out that propagation through 
such a medium exhibits many new and interesting properties, such as power oscillations and 
^ ■ birefringence. The equation of propagation takes the form 



where tp(x,z) represents the envelope function of the amplitude of the electric field, z is a 
scaled propagation distance, and V(x) is the optical potential, proportional to the variation 
in the refractive index of the material through which the wave is passing. A complex V 
corresponds to a complex refractive index, whose imaginary part represents either loss or 
gain. In principle the loss and gain regions can be carefully configured so that V is PT 
symmetric, that is V*(x) = V(—x). There is also a non-linear version of this equation, 
arising from sufficiently intense beams, where there is an additional term proportional to 
|V'| 2 V'- However, for the purposes of this paper we shall limit ourselves to the linear case. 
A model system exemplifying some of the novel features of beam propagation in PT- 
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symmetric optical lattices uses the sinusoidal potential 

V = Vq [cos(27rx/a) + iX sin(27rx/a)] 

This model has been studied numerically and theoretically in Refs. [9], Qjl, 13 1. The prop- 
agation in z of the amplitude if)(x, z) is governed by the analogue Schrodinger equation!], 
which for an eigenstate of H, with eigenvalue (3 and ^-dependence if) oc e~ l ^ z reduces to the 
eigenvalue equation 

- if;" - V [cos(27rx/a) + iA sin(2vrx/a)] ip = f3if; . (2) 

It turns out that these eigenvalues are real for A < 1, which corresponds to unbroken PT 
symmetry, where the eigenf unctions respect the (anti-linear) symmetry of the Hamiltonian. 
Above A = 1 complex eigenvalues begin to appear, and indeed above A ~ 1.77687 all 
the eigenvalues are complex[l~6T]. Clearly one would expect oscillatory behaviour of the 
amplitude below the threshold at A = 1 and exponential behaviour above the threshold, but 
the precise form of the evolution at A = 1 is less obvious. At first sight one would expect 
linear growth because of the appearance of Jordan blocks associated with the degenerate 



eigenvalues that merge at that value of A, but, as Longhi[l2j has emphasized, this behaviour 



can be significantly modified depending on the nature of the initial wave packet. 



In a previous paper[17| we approached this problem by explicitly constructing the Bloch 
wave-functions and the associated Jordan functions corresponding to the degenerate eigen- 
values and then using the method of stationary states to construct the dependence. We 
found that the explicit linear dependence arising from the Jordan associated functions is 
indeed cancelled by the combined contributions from the non-degenerate wave-functions 
and were able to understand how this cancellation came about. In the present paper we 
approach the problem from a different point of view by revisiting the complementary pertur- 
bative calculation of Longhi[12|. In Section 2 we briefly recapitulate how the spectrum and 
eigenfunctions are calculated. Then in Section 3, which forms the main body of the paper, 
we give an explicit expression for the first-order contribution and carry out the second-order 
calculation in detail. This enables us to investigate the saturation phenomenon for a variety 
of different inputs. Finally, in Section 4 we give a brief discussion of our results. 



II. BAND STRUCTURE AT THRESHOLD 

At the threshold A = 1, the potential V in Eq. 02]) becomes the complex exponential 
V = Vq exp(2i7TX / a) , for which the Schrodinger equation is 

- ip" - V exp(2iirx/a)i/> = pi/). (3) 

This is a form of the Bessel equation, as is seen by the substitution y = y Q exp(i7rx/a), where 
yo = (a/7r)v%, giving 

«"0 + v|-to» + ^ = O, (4) 

where q 2 = P(a/n) 2 . Thus the spectrum is that of a free massive particle, shown in the 
reduced zone scheme in Fig. 1, and for q = ka/ir not an integer the solutions ipk{x) = Iq{y) 
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and ip_ k (x) = I- q (y) are linearly independent, and have exactly the correct periodicity, 
ip k (x + a ) — e lha ijjk {%) , to be the Bloch wave-functions. It is important to note, however, 
that because the original potential is PT-symmetric rather than Hermitian, these functions 
are not orthogonal in the usual sense, but rather with respect to the PT inner product, 
namely 

J dxi)_ k {x)^ k >{x) = S kk > J dxi)„ k {x)4) k {x), (5) 
However, for q = n, a non-zero integer, I n (y) and I- n {y) are no longer independent. In 




FIG. 1: Band structure for A = 1 in the reduced zone scheme. The Bloch momentum k is plotted 
in units of Tt/a and the eigengvalue j3 in units of (a/-7r) 2 . 



that case the Bloch eigenfunctions do not form a complete set, and we must search for 
other functions, still with the same periodicity, to supplement them. These are the Jordan 
associated functions, which we denote by (p k (x) = Xn(y)- They may be defined as derivatives 
of the eigenfunctions with respect to (3, and satisfy the generalized eigenvalue equation 



d 2 d 

y -n + y-r - \y + n ) 

dy dy 



Xn{y) = l n (y), 



(6) 



The crucial feature of the Jordan functions is that because of this latter equation they 
naturally give rise to linear growth in z, provided that they are excited: 



e~ iHz ip r 



e~ il3rZ ((p r - izip r ). 



(7) 



However, as was found numerically in Ref. [12j, and explored further in Ref. [17], this 
natural linear growth may become saturated due to the contributions of neighbouring Bloch 
functions, which are closely correlated with those of the Jordan functions. 
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III. PERTURBATION THEORY 

The analysis of Ref. [l7[ approached the problem from one point of view, in which the 
interplay between the contributions of the Bloch eigenfunctions and the Jordan associated 
functions was made explicit. A complementary way of looking at things, which does not 
separate these two contributions, is to use the perturbative expansion, which instead em- 
phasizes the contributions of the free propagation and the corrections brought about by 
the potential. The general framework for an expansion of ip{x, z) in powers of Vq, namely 
ip(x,z) = Y^°=o KTVv^j z )i has been given in Ref. 12), along with an approximate form of 
the first-order term ipi(x, z) for the case q$ = —1 and w large. In this section we generalize 
this calculation by obtaining analytic expressions for both the first- and second-order terms 
for general q and w. Of course, this can only be used as a guide because there is no guaran- 
tee that the expansion converges, nor that the large- z behaviour of the complete amplitude 
can be extracted from the behaviour of the truncated series. We will take as our input a 
Gaussian profile of the form 

4>(x,0) = f(x) = e-C*/*) 2 -*** (8) 

with offset fco and width w. The zeroth-order term, ipo(x,z), is just the freely-propagating 
wave-packet 

J.( T -\ _ W ik (x-k„z) -(x-2k Q z) 2 /(w 2 +4iz) /q\ 

W ' ] ~ j(w 2 + Aiz) ' [ ) 



while the first-order term, ipi(x, z), is given by 

ijji{x,z) = -i J dkf{k)e i{k+ks)x 

where 



dy e~ lK y e 



ik 2 y i(k+k B ) 2 (y-z) 



(10) 



/(*) = -^e'^ 2 ^ (11) 

is the Fourier transform of f(x) of Eq. (jSJ) and = 2ir/a is the width of the first Brillouin 
zone. We can reverse the order of integration in the expression for ipx, performing the 
(Gaussian) k integration first, to obtain 

ibAx.z) = -l „ W e \ik B x-\iklz~\{k + lk B f W 2 x (12) 

v ' y/{w 2 + 4iz) v ' 

dy e -{ 2k By-(x+kBz)+\iw 2 {k Q + \k B )) 2 /(w 2 +4iz) 

The y integration is then also a Gaussian integration over a finite range, giving the result 



„• p \ik B X- \ik 2 B Z- \ (fco + \ kg ) 2 W 2 

where 



ipx(x,z) = -i^— e -2^BX—^k B z—^+- 2 k B) w (erf^) + er% 2 )) , (13) 
Aks 



k B z + x - \iw 2 (ko + \k B ) 
Vl = v/(w 2 + Aiz) 
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n-2 



k B z — x + \iw 2 (kQ + \k B ) 



^{w 2 + Aiz) 

For the purposes of considering large w, it is convenient to rewrite these in the form 

2koZ 



(14) 



x 



Vi 

m 



^(w 2 + Aiz) 
2z(k B + fa) 



- -i(k + ^k B )y/(w 2 + Aiz) 



^ + Aiz)- + l i{ko + l kBMw2 + 4iz) 



(15) 



The case kn 



to = —^^B, i.e. go = — 1, is clearly very special, since in this case the second 
terms in the contributions to rji and r]2 vanish, so that we get the simple expressions r/i = 
(x + k B z)/ y/(w 2 + Aiz) and 7/ 2 = —(x — k B z)/ \/{w 2 + Aiz). In that case, as long as w 2 ^> Az 
the arguments may be treated as effectively real, and each error function behaves like a sign 
function of its argument (see Fig. 2(a)), so that the sum of the two behaves like the step 
function 6(k B z— \x\). This is the function Q(x/(k B z)) of Ref. T3 |. In this case the qualitative 
features of the perturbative calculation are in complete agreement with the spreading of the 
wave-function in Fig. 3 of that paper, and the saturation 1 of VWx- However, in this treatment 
there is of course no mention of whether or not any Jordan functions are excited. 
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FIG. 2: Characteristic behaviour of error functions whose arguments are (a) real or (b) have a 
large imaginary part. In (a) we plot erf(4 + y/w) + erf(4 — y/w) with w = 80. In (b) we plot 
w e~ w | erf (4 + y/w — iw) + erf(4 — y/w + iw)\ 



The same expressions in Eqs. (1T3"]) and (TT51) can also be used for the cases q = and 
go = 1 in the limit of large w. In each case the arguments rji and r]2 now have a large 
imaginary part. In that situation the modulus of the erf has a narrow peak where the real 
part vanishes (see Fig. 2(b)). Thus the result consists of two narrow rays, which are centered 
on x = 2koz and x = 2(k B + ko)z. Here we have the seeds of the birefringence first observed 
in Ref. For the case ko = qo = 0, the two rays are centered on x = and x = 2k B z, 
while for the case go = lj or = the two rays are centered on x = k B z and x = 3k B z. 

We now go on to second-order perturbation theory to investigate the behaviour of ^(^ z), 
which is given by[l2T| 

fo(x,z) = - [ dkf(k)e i( - k+2kB) f drj [ d£ e -^ 2 z-4ik b (k+k B ) v +ik B (k B +m (i 6 ) 
J Jo Jo 



4>i{x, z) grows initially like z for small z. 
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Again the k integration is a Gaussian, which leaves finite-range Gaussian integrations over 
£ and rj. It is convenient to change the integration variable £ to y = — £. The integration 
over 1/ then yields the expression 



$2(3,*) 



f ^^(erf(a) - erf(&)) e - 2ifc ^ e § ^-i*I*-X a2 , 
Jo 2 



where we have written k = — {\k B + A), and a and 6 are given by 

b = 



4k B r) — x — k B z — ^iw 2 A 



^(w 2 + Aiz) 
3k B (2r) -z)-x- \iw 2 A 
yj(w 2 + 4iz) 

The final 77 integrations are then of the form 



(18) 



I v = J drj eri(cir] + C2)e C3V 

e C3V erf{ Cl r] + c 2 ) - e c3( C 3-4 ClC2 )/(4 C 2) erf + ^ _ C3 /( 2ci ))] . (19) 



1 

c 3 



Thus in principle ?/> 2 is expressible in terms of eight error functions. However, it turns out 
in practice that only six are involved. 

In the case k = — /c.b/2, or q = — 1, when A = 0, the arguments of the error functions 
are such as to give a plateau in \ip 2 \ between x = Sksz and x = —ksz, i.e. a widening of 
the beam to the left of ipo, and a much smaller peak, centered around x = 3k B z, that is, 
a second weak beam to the right. More importantly, the second-order contribution again 
shows no sign of the linear growth 2 in z naively expected from the excitation of Jordan 
associated functions. 

In the other case, q = 0, corresponding to A = — there are three peaks, centered 
around x = 0, x = —2k B z and x = 4k B z, representing a further splitting of the initial beam. 
Both cases are illustrated in Figure 3. 




FIG. 3: \ip2 (•£) z)\ versus x for z = 50. (a) qo = —1, (b) qo = 0. The parameters are: a = 1, Vq = 2 
and w = 6tt. 



2 In fact tp2(x, z) is proportional to z 2 for small z. 
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IV. DISCUSSION 

Of course one needs to treat the results of perturbation theory with caution. In gen- 
eral terms we have no proof that the perturbation series converges, and in particular the 
asymptotic behaviour in z of a few terms of the series does not necessarily give the correct 
asymptotic behaviour of the entire sum. Nonetheless this series does appear to give reliable 
results. For the parameters used by Longhi in Ref. [l2|, with Vo = 0.2, first-order perturba- 
tion theory already reproduces the numerical results very well, and the second-order term 
gives only an extremely small correction. 

In all of our calculations we have not allowed z to become too large, restricting it by the 
condition z <C w 2 /4, in which case saturation is an inbuilt feature of perturbation theory. 
In fact it was shown in Ref. [l7[ using the method of stationary states that if one goes to 
much larger values of z the amplitude so calculated begins to grow again but this ultimate 
resumption of linear growth is not physical, because it corresponds to the situation where 
the beam has widened beyond lateral limits of the optical lattice. 

It is an elegant feature of the perturbative expansion that the different types of possible 
behaviour of the beam - spreading or splitting into two or more beams - arise from very 
simple properties of the error function depending crucially on the offset k$. In the first case 
the arguments are essentially real, and the error functions behave like sign functions, while 
in the second case there is a large imaginary part, and the moduli of the error functions 
behave instead like narrow peaks. 
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